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We model buoyancy-driven convection with chimneys - channels of zero solid fraction - in 
a mushy layer formed during directional solidification of a binary alloy in two dimensions. 
A large suite of numerical simulations is combined with scaling analysis in order to 
study the parametric dependence of the flow. Stability boundaries are calculated for 
states of finite-amplitude convection with chimneys, which for a narrow domain can be 
interpreted in terms of a modified Rayleigh number criterion based on the domain width 
and mushy-layer permeability. For solidification in a wide domain with multiple chimneys, 
it has previously been hypothesized that the chimney spacing will adjust to optimize the 
rate of removal of potential energy from the system. For a wide variety of initial liquid 
concentration conditions, we consider the detailed flow structure in this optimal state and 
derive scaling laws for how the flow evolves as the strength of convection increases. For 
moderate mushy-layer Rayleigh numbers these flow properties support a solute flux that 
increases linearly with Rayleigh number. This behaviour does not persist indefinitely, 
however, with porosity-dependent flow saturation resulting in sub-linear growth of the 
solute flux for sufficiently large Rayleigh numbers. Finally, we consider the influence of 
the porosity dependence of permeability, with a cubic function and a Carman-Kozeny 
permeability yielding qualitatively similar system dynamics and flow profiles for the 
optimal states. 



1. Introduction 

Solidifying binary alloys occur in a wide variety of geophysical, geological and industrial 
settings. In many contexts a mushy layer is formed as a result of morphological instability 
of an advancing solid-liquid interface. A mushy layer consists of a reactive porous medium of 
solid dendrite s in local thermodynamic equilibrium with the interstitial concentrated fluid (e.g., 
Worster 2000). An important geop hysical example is the growth of sea ice via the solidification 
of salty water (|Feltham et al\\200$ ), with ~ 10 6 -10 7 squa re kilometres of the Arctic and Antarc- 
tic oceans freezing and melting each year (Comiso 2010). The drainage of b rine from growing 
sea i ce modifies the thermal and mechanical properties of the ice (see e.g., iPetrich fc Eickea 
2010;, for a recent review) and plays a significant role in the formation of dense w ater masses 
that contribute to the ocean thermohaline circulation (see e.g.. iBrandon et a/.ll2010l . for a recent 
review). I n industrial setting s, defects in metal castings can arise from compositional hetero- 
geneities (jCoplev et a/.ll 1970) 1, and so it is important to understand the evolution of solidifying 
alloys. Here we focus on the influence of buoyancy-driven convection on the drainage of solute 
from a growing two-dimensional mushy layer. 
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Buoyancy-driven convection of th e interstitial flu id in a mushy layer occurs when the density 
gradient is convectively unstable (sec Worstcr 1997, for a review). The convective flow transports 
interstitial fluid into regions of differing concentration, which results in either local dissolution 
or growth of the solid matrix in order to maintain the composition in a state of local thermody- 
namic equilibrium. In regions where dissolution decreases the solid fraction, the permeability is 
increased leading to flow focussing. The nonlinear development of this convective-flow instabil- 
ity forms drainage channels of zero solid fraction, or chimneys, through which buoyant plumes 
drain from the mushy layer into the neighbouring liquid. Experimental observations are consis- 
tent with chimney formation occurring for sufficiently large values of an appropriate mushy-layer 
Rayleigh number, which characterizes the strength of buoyancy compared to viscous and ther- 
mal dissipation. For directional solidification, where a material sample is translated between 
hot and cold heat exchangers, chimneys form for sufficiently small solidification rates V so that 
the resulting steady-st ate mushy layer thi c kness h is large enough to generate a sup ercritical 
Rayleigh number (e.g., iPeppin et atll200Sl ; IWhiteoak et aZj|2008T ). IZhong et all (|2012f ) showed 
that the horizontal wavelength provides a further constraint on chimney formation in cells of 
finite width. Chimney formation also occurs after the mushy-layer thi ckness h exceeds a crit ical 
value during transient solidification from a fixed cold boundary f see lAussillous et a/.ll2006l . for 
a review). Transient growth experiments also indicate that the pattern of chimneys continues 
to evolve over time, with the extinction of flow in certain chimneys le ading to an increase in 
the mean spacing betwee n chimneys as the mushy layer grows thicker (Wettlau fer et a/"1 ll997: 
ISolomon fc Hartley! Il998f ). 



Theoretical studies have identified critical Rayleigh numbers for the onset of convection under 
a w ide range of growth conditions by ap plying linear and weakly-nonlinear stability analyses 
fsee lWorster!ll997l : lGuba fc Worstedl2010l . for a review and a more recent summary) . The forma- 
tion of chimneys presents additional challenges for modelling, which must describe and resolve 
the change in flow dynamics between the interior of the porous mushy layer and the purely 
liquid region in the chi mney. Single-domain methods, such as the enthalpy method approach of 
Kat z fc Worsted (|2008T ). have simulated solidification with small numbers of chimney- like fea- 
tures in a vari ety of configurat i ons (s ee lZnong et aHl20I2l . for a summary). In particular, the 
simulations of lKatz fc Wo rstcr (2003) f° r directional solidification in a quasi-two-dimensional 
Hele-Shaw cell suggest that the mean chimney spacing scales in proportion to the depth of the 
mushy layer. An alternative approach takes advantage of the fact that chimneys are narrow 
compared to th eir depth thereb y allowing the use of lubrication theory to parametrize the flow 
in the chimney. IWorster ] (|1991 ) considered axisymmetr ic flow with an isolated chimney, whilst 



ISchulze fc Worsted (|199gl ) and I Chung fc Worsted (|2002T ) developed a model for two-dimensional 



steady flow with a periodic array of chimneys with imposed spacing. Theories that require the 
unknown chimney spacing to be imposed in advan ce of the calculatio n have several advantages 
and varying degrees of sophistication. For example, IWells et al\ (|2010h combined this constraint 
with a hypothesis of optimal potential energy fluxes, so that the system selects the chimney 
spacing that yields the maximal rate of drainage of potential energy to efficiently drive the 
system towards thermodynamic equilibrium. This hypothesis predicts chimney spacings pro- 
portional to mushy layer depth and a solute flux that increases linearly with Rayleigh number 
over the simulated range. The resulting solute flux scaling predicts initi al desalinization rates 
during sea ice growth that are consistent with experimental observations l] Wells et ffl/.ll20liT ). An 
imposed chimney spacing approach also allows the development of an asymptotically-reduced 
model that reveals the flow structure for mushy layers of constant perm eability and assumed 
constant background temperature gradient (|Rees Jones fc Worstedl2013l ). 



We here explore in detail the framework of IWells et al\ (|2010T ) to determine the leading order 
structure of the flow and its parametric dependence for differing system properties. We consider 
convection during the directional solidification of a two-dimensional mushy layer with a periodic 
array of chimneys. In fj2] we describe the theoretical model and numerical methods employed. 
In >J3] a broad array of numerical simulations are used to demonstrate the evolution of the 
system dynamics and stability as the chimney spacing and vigour of convection change. We 
then investigate the dynamics of states with optimal potential energy flux in Sj4] and for a 
variety of liquid concentrations determine asymptotic scaling laws for flow properties as the 
strength of convection increases. We consider the influence of the permeability variation with 
solid fraction in section Sj5] before concluding in $6] 
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Figure 1. Schematic illustration of the model configuration for directional solidification. Liquid 
of far-field salinity Co and temperature Too is translated at velocity V between hot and cold heat 
exchangers, and partially solidifies to form a mushy region of thickness h(x, t). A eutectic solid 
forms at the lower heat exchanger with temperature Te and concentration Ce at the mush-solid 
interface. In this geometry, with gravitational acceleration — <?k aligned perpendicular to the 
heat exchangers, solidification releases buoyant fluid of relative density p g/3(C — Ce) to drive 
the convective flow, where p is a constant reference density and /3 the constant haline coefficient. 
We consider a periodic array of chimneys of dimensional width 2a(z, t) and imposed spacing I, 
with the dashed outline indicating the simulation region. Buoyant solute-depleted plumes exit 
the mushy layer via chimneys, and there is a weak return flow in the overlying liquid as indicated 
by arrows in the liquid region. Along with other systems, this geometry is consistent with the 
solidification of a trans-eutectic aqueous ammonium chloride solution cooled from below, where 
cooling at the base of the mushy layer lowers the liquidus temperature and results in lower solute 
concentration of the interstitial fluid. 



2. A model of non-linear convection with chimneys 

Building on the previous models of iSchulze fc Worster! (|l998l ) and lChung fc Worsted IpOOl ). 
we apply the framework of IWells et al.\ (|201oh to describe two-dimensional directional solid- 
ification, as illustrated in figure [T] A two-component mixture of liquid concentration C and 
temperature T is translated at a velocity V between hot and cold heat exchangers, and forms 
a mushy region of thickness h(x,t). In the illustrated setting, solidification within the mushy 
region lowers the local concentration of the liquid, and hence provides the buoyancy to drive 
convection. This would occur in aqueous NH4CI solidified from below and the dynamics and 
thermodynamics are ostensibly the same as ice forming above a region of salt water. We inves- 
tigate the behaviour of a periodic array of chimneys within this system. For a given chimney 
spacing /, we calculate the resulting solute flux. 

The thermodynamic evolution of the mush y region is treated within the framewo rk of so- 
called ideal mushy layer theory (for example, sec Worster 2000; Schulzc & Worster 2005, for a full 
discussion of the underlying physics). It is assumed that the specific heat capacity c v and thermal 
diffusivity k are constant across both solid an d liquid phases, the li quid has kinematic viscosity 
v, and L is the latent heat of fusion. Following [Schulze fc Worsted (|l998f ). we scale velocities by 
V, and length and timescales via the thermal diffusion scales k/V and k/V 2 respectively. The 
temperature and concentration within the mushy layer are constrained by the condition of local 
thermodynamic equilibrium, with the liquid solute concentration maintaining a state of local 
phase equilibrium by either dissolving or crystallizing the solid matrix faster than any timescale 
of macroscopic transport. Hence the liquid concentration is slaved to the temperature via the 
liquidus curve T = Tl(C) = Te + T(C — Ce), where F is constant. Both temperature and liquid 
concentration can then be characterised by the local dimensionless temperature 

g _ T - T L (C ) _ C -Co ^ ^ 

Tl{Co) — Te Co — Ce 

where Co is the concentration in the bulk liquid layer (see Fig.[T}. The dimensionless temperature 
6, solid fraction <j>, Darcy velocity u and pressure p then satisfy 



u = 



-RmII(Vp + 6>k) , 



V ■ u = 0, 



(2.2a,6) 
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J^- J^) [(1-^)0 + C0]+u-V0 = 0, (2.4) 

representing conservation of momentum, mass, energy and solute. We have assumed that the 
fluid density depends linearly on concentration, that solute diffusion is negligible compared to 
thermal transport, and that the mushy layer has a porosity-dependent dimensional permeability 
IT = IloIl((f>), where ITo is a reference permeability. We will consider two different functional 
forms of the permeability. A cubic permeability function 

U(<j>) = (1 - 0) 3 (2.5) 

is consistent with earlier studies (|Schulze fc Worsted 1 19981 ; IChung fc Worstcr 2002) and bears 
close rese mblance to the permea bility dep endence observed for se a ice over a range of solid 
fractions (|Petrich fc Eickerj2010l ). Following iKatz fc Worsted" feOOS) we also consider a selection 
of simulations in fusing a modified Carman- Kozeny permeability 



i 



(2.6) 



This approximates the usual Carman-Kozeny form for larger (j> and the factor e = 12ilo/d 2 
accounts for sidewall drag in narrow Hele-Shaw cells of width d and removes a singularity as 
<j> — > 0. There are six dimensionless parameters governing the system, 

_ gP(C - C E )n L Cs-C 

Rm - uV ' S= c p F(Co-C E y C= C^C- E ' (2 ' 7) 

Da -^' 9 °°~ T E -T L {c y (2 - 8) 

where the mushy layer Rayleigh number Rm describes the ratio of buoyancy to dissipation, 
and the Darcy number Da characterizes the mushy layer permeability. The Stefan number 
S, concentration ratio C, in which Cs is the solid concentration, and scaled temperature 8oo 
characterize the imposed thermodynamic conditions. The dimensionless chimney half-spacing is 
A. Taking the curl of (|2.2b ) yields the vorticity equation 



Vtf) • VIA 



where the streamfunction tp is defined by u = (—dip/dz,dip/dx), thereby satisfying (\2.2b ). 

Flow features in the chimneys and overlying liquid can occur on smaller length scales than 
features in the mushy region, and so a considerable increase in computational effort w ould 
be required to fully resolve the liquid regions numerically (e.g. IChung fc W orstcr 2002|). In 
this study, we focus on dynamics dominated by the mushy region, and utilize asymptotic and 
perturbation approximations to describe the influence of the overlying liquid and chimney flow 
via boundary conditions for the mushy region. The boundary conditions at the mush-liquid 
interface z — h describe coupling to the overlying liquid. If solute diffusion is neglected, the 
liquid region has uniform concentrat ion Co outside of the compositional plumes exiting the 
mushy layer (ISchulze fc W orstcr 1991), and we assume there is constant pressure at the mush- 
liquid interface (|Emms fc Fowled 1 1994T ). The mush-liquid interface advances until it eliminates 
constitutional supercooling, so that the free-boundary position z — h(x,t) is determined from 
the condition of marginal equilibrium T = Tl(Co). Combining continuity of temperature and 
solute concentration, the above assumptions yield 

6» = 0, <j} = 0, n-Vtp = 0, at z = h, (2.10a,b,c) 

where Darcy's law (|2.2b ) has been used to write the pressure condition in terms of ip. The final 
boundary condition at the mush-liquid interface comes from continuity of normal heat fluxes 
n ■ VT|^ = 0. Balancing advection and diffusion of heat across a thermal boundary layer where 
isotherms have curvature V • n yields 

n-V0 = 0oo[V-n-(u-k)-n], {z = h). (2.11) 

(see IChung fc W orstcr 2002L for a full derivation.) 
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We describe the flow in the chimney by combi ning features of previous analyses (jSchulze fc Worsted 
1998; Chung fc Worstcr 2002; Wells et al. 2010), taking advantage of the observation that chim- 
neys are narrow with a ft to derive representative singular-interface conditions that are 
applied at x — 0. Note that the singular-interface approximation can formally be justified via a 
stretched co-ordinate transform of the form x = X(x — a)/(A — a) in the horizontal, and gives rise 
to corrections of 0{a) <C 1 in the governing partial differential equations. Neglecting these small 
corrections of O(a) recovers the structure of the original partial differential equations with x re- 
placing x, but with the mushy layer occupying < x < A. However, in what follows we will drop 
the annotation and proceed with notation where the mushy region occupies < x < A with the 
chimney described by boundary conditions applied at a singular interface at x = 0. The chimney 
width a(z, t) then acts as a parameter in the boundary conditions al ong the singular interfac e 
at x = 0. Fluid flow in the chimney is described by lubrication theory (|Chung fc Worster]|2002f ). 
yielding 

g + (-0), <»..»> 

where a quadra tic Polhausen approximat ion has been used to determine the pre-factor in the 
buoyancy force (Schulzc & Worstcr 1998). The mass-flux boundary condition (|2.12[) is applied 
to the vorticity equation Q2.90 . The term adijj/dx in (|2.12|l results from assuming that, in the 
limit of vanishing solid fraction, the velocity is continuous at the mush-l iquid interface that 
forms the chimney wall (|Chung fc W orstcr 2002T). [Le Bars fc Worstcr (200(1) demonstrated that 
an improved porous medium boundary condition imposes continuity of u at a small distance 
S of order the pore scale inside the porous medium. However, for simplicity we here retain the 
boundary condition (|2.12|l which can be viewed as the leading term in a Taylor expansion for 
5< 1. 

The heat conducted into the chimney wall balances the heat flux advected along the chimney, 
yielding 

l = *t <"»)■ ^ 

The heat flux boundary condition (|2 . 1 3|) is applied to the heat equation (|2.3[) . As demonstrated in 
appendix(X] this boundary condition relies on the approximation atp/h <C 1. This approximation 
was previously justified in the asymptotic limit Rm 4 ^ 3 Da <C 1 bv lSchulze fc Worsted l|l998f ). 
The simulations presented here extend over a wider range of parameter space, and so we proceed 
under the ansatz aip/h <C 1 and treat (|2.13|) as a leading order perturbation expansion. Our 
subsequent simulations reveal a/h < 0.06 and tp < 1, so this is a self-consistent assumption. 

The chimney width a(z, t) is determined from marginal equilibrium, which for a free boundary 
with net outflow, yields 

riO rift 

f-g+u-Ve = 0, (x = 0), (2.14) 

ijSchulze fc Worsted l2005h . The chimney width a(z,t) adjusts in order to satisfy the condi- 
ti on (2.14). This is t he tim e-dependent generalization of the marginal equilibrium condition 
of lChung fc Worsted (|2002l ). 

The remaining boundary conditions on the lower and right-hand boundaries of the mushy 
region are 

9 = -l, ip = 0, at 2 = 0, (2.15) 

d8 

— = 0, ip = 0, at x = A. (2.16) 

ox 

corresponding to no mass flux across the lower heat exchanger which is fixed at the eutectic 
temperature, and symmetry conditions at x = A. 

We pause here to consider some restrictions regarding the current modelling framework using 
ideal mushy layer theory. The use of Darcy's law to describe porous medium flow assumes a 
low Reynolds number Re p = US/u = O(l) within the pore space, where U is an approximate 
dimensional velocity scale and 5 an approximate lengthscale of the pores. In addition, the con- 
tinuum conservation equations are justified by volume averaging across solid and liquid phases, 
and assume that all macroscopic flow features within the mushy layer occur on lengthscales 
larger than the pore scale. We return to these conditions in £14.21 
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2.1. Numerical methods 

We solved the system (|2.3p - (|2.16[) numerically using second-order finite differences for the time- 
dependent problem, complemented by an arc-length continuation method that traces both stable 
and unstable solution branches in order to verify the bifurcation structure. A brief summary 
follows, with a more detailed elaboration in Appendix [Bl 

Semi-implicit Crank-Nicholson time-stepping was used to update the heat and concentration 
equations (|2.3[) and (|2.4|l . while the vorticity equation (|2.9[) yields a nonlinear Poisson equation 
for ip. After time discretization, the spati al derivatives of 1)2. 3[1 and f[2T9l) yield elliptic systems 
that were solved using multigrid iteration ((Adams 1989; Bri ggs et a/.ll200(J ). The free boundaries 
a(z, t) and h(x, t) were updated via relaxation, with a(z, t) updated at each value of z to reduce 
the error in (|2.14|l . A corresponding free-boundary problem for h(x,t) using (|2.10b ) yielded 
an unstable sche me when used in conjunc tion with the boundary layer approximation (|2.11[) . 
Hence, following ISchulze fc Worsted (| 19981 ). a one-parameter shape 

il>c fl — cosh u (A — x)] ., 

is enforced, which has a thermal boundary layer in order to remove a temperature singularity at 
the chimn ey top. We use a th e rmal- boundary-layer width 1/fx — Rm~ 2//3 following the scaling 
analysis of Schulzc fc Worsterl (|1998l ). Note that (|2.17|l is an approximation to h(x,t) and does 
not provide an exact solution of the free boundary problem. We choose to apply the boundary 
conditions so that the condition of marginal equilibrium (|2.10b ) is enforced exactly at all times 
by setting 9 — Oaiz — hasa boundary condition on (|2.3[) . and then ho is updated to minimize 
the least-square residual in the heat flux condition (|2.11|l . In this manner, errors from the two 
approximations (|2.11|l and (|2.17[) are incorporated in one condition. The time-dependent initial 
value problem was integrated to a steady state with initial co nditions given either by an analytic 
solution for solidification with no fluid flow (e.g. lWo rstcr 199l]) or by continuation from a previous 
solution with different parameters. In addition to other flow quantities, we diagnose the solute 
flux from the mushy layer into the fluid per unit width. The combination of conservation of 
solut e with continuity of C at a freezing interface with outflow results in cj> = at the chimney 
wall (jSchulze fc Worsterll2005l ). and further noting that — </> = at z = h, the resulting solute 
flux reduces to a line integral over the chimney boundary, with absolute magnitude 



F a = i / 9{u-k)ndl. (2.18) 

/{»=a00} 



A 



In steady state, integrating ()2.4|l over the area of the mushy region, using the divergence theorem 
and applying the boundary conditions (|2.10[) . (|2.15[) and (|2.16[) yields the alternative expression 

F s = ±J [(l-£)0 + C4>U o dx > ( 2 - 19 ) 

which is adopted here for numerical convenience. 

To complement the time-dependent solutions, an arc-length continuation scheme was applied 
to tr ace both stable and unstable steady states, and verify the bifurcation structure (I Keller! 
Il977f l. To reduce the computational cost associated with the size of the Jacobian matrices 
required for arc-length continuation, we project the finite-difference solutions of (|2.3[1 - (|2.16[1 
onto a basis of Chebyshev polynomials yielding a low-order pseudo-spectral representation for 
use in the predictor-corrector scheme. The resulting updated Chebyshev projections were then 
used as initial conditions in an iterative time-independent version of the finite difference code 
for post processing (see Appendix iBl for further details). 



3. Flow dependence on chimney spacing 

We begin by considering the ev olution of flow dyn amics as the chimney spacing is varied, illus- 
trating key features described by I Wells et al\ lj2010l ) to provide context for the later discussion. 
To demonstrate the system dynamics, figure [2] shows a characteristic example of the variation 
of the dimensionless solute flux per unit width F 3 with chimney half-spacing A. The system 
has two stable states showing hysteresis: a lower branch of no flow for A < X u (solid line with 
triangle symbols) and an upper branch of chimney convection for A > A c (solid black curve). 
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Figure 2. Calculated steady-state solute flux per unit width F s as a function of chimney 
half-spacing A. The steady states show hysteresis, with a stable upper branch of convection with 
chimneys for A > A c (black curve) and a lower branch of no flow for A < \ u (blue line with 
triangle symbols) separated by an unstable solution branch (red dashed curve) . A fold bifurcation 
is observed at (A c , F c ), and the optimal flux F is achieved with chimney half-spacing A . These 
fluxes are calculated for Rm = 95, C = 5, S = 5, 9oo = 0.4 and Da = 5x 10 -3 . 



The arc-length continuation method also reveals an intermediate unstable branch of chimney 
convection (dashed curve), with the two states of chimney convection becoming extinct via a 
bifurcation at (X C ,F C ) as A is reduced. For large values of Rm this corresponds to a saddle- 
node bifurcation. For much smaller Rm, close to the linear critical value of Rayleigh number, 
the arc-length continuation suggests a singular Jacobian matrix at the turning point consistent 
with a more complex bifurcation structure involving flow without steady-state chimneys (see 
for example the bifurcation to a state of convection without a chimney shown in figure 2b of 
Wells et al. 20l3). However, we emphasize that, for all simulations considered, the solute flux 
shows the same qualitative turning point structure as illustrated in figure[2] with steady states of 
chimney convection ceasing to exist for A < A c . Failure of convergence of the arc- length scheme 
prevents continuation of the unstable branch to lar ger A. The s i mulat ions considered here show 
that this form of hysteresis behaviour identified by I Wells et al\ (|2010f ) persists for a wide range 
of concentration ratios C and Rayleigh numbers Rm. 

The solute fluxes illustrated in figure [2] have an optimal value at A = A . Figure [3] shows the 
corresponding profiles of temperature (colour scale in panel a), solid fraction (colour scale in 
panel b) and chimney width (panel c) in this optimal state, along with streamlines of Darcy 
velocity (magenta curves in panel a) and net fluid flux q = u — k relative to the heat exchangers 
(magenta curves in panel 6), wh ere k is a uni t vecto r in the ^-direction. These profiles exhibit the 
same flow structure observed bv lWells et al\ ()2010l ). with inflow at the top boundary, circulation 
in an order-one aspect ratio convective cell, and then outflow into the chimney at the left side of 
the domain. At the large Rayleigh number of this simulation we observe a tapered chimney shape 
which is narrow at th e base, differing slightly from the steep-sided chimneys observed at low 
Rayleigh number (see IChung fc Wo rstcr 2002], and also confirmed using our current method). 
In ^4] we determine a transition in scaling regimes for the chimney between dominance of plug 
flow and Poiseuille flow, which may underlie this difference in shape. 

We now examine in detail the flow dynamics as the Rayleigh number and concentration ratio 
are varied. For fixed Rayleigh number, the bifurcation at (X C ,F C ) defines a stability boundary 
for the existence of the observed chimney-convection state in cells of confined width. Figure [4] 
traces the stability boundaries A c (Rm) as a function of Rayleigh number for a wide range of 
concentration ratios with C > 2. Chimney convection is stable above each of the curve sections 
plotted in the (A, Rm) plane in figure \4ja. All concentration ratios used are consistent with the 
same underlying pattern, with chimney convection suppressed for small Rayleigh numbers o r 
small chimney spacings. This behaviour is consistent with the experiments of lZhong et al\ (|2012f ). 
Note that convergence issues in the arc-length continuation method prevented us from extending 
the range of these curves: the failure of convergence of the corrector scheme is responsible for 
failure at larger Rayleigh numbers, and for the smaller Rayleigh numbers, due to the existence 
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Figure 3. Steady-state mushy layer profiles with an optimal chimney half-spacing A = 0.194, 
and all other parameters identical to figure [2] In panel (a) the colour scale shows dimensionless 
temperature 9{x,z) and solid black curves show isotherms, whilst in panel (b) the colour scale 
shows solid fraction <j>(x,z), with black solid contours of constant (j). Magenta curves show 
streamlines of Darcy velocity u in panel (a) and net fluid flux q — u — k in panel (6), with 
inflow at the upper boundary of the mush, and outflow into the chimney at the left boundary. 
Panel (c) illustrates variation of the chimney width a(z) with height. 




Figure 4. (a) Stability boundaries A c (Rm) for a wide variety of concentration ratios C as given 
in the legend. For a given value of C, stable states of chimney-convection exist everywhere above 
the plotted section of curve (i.e. A > A c (Rm)) in the (A,Rm) plane, (b) Scaling of stability 
boundaries illustrated by re-plotting 1/Rm versus A c for all data in panel (a). Line styles and 
symbols follow the legend in panel (a). In both plots, all other parameters S — 5, 8^ = 0.4 and 
Da = 5x 10~ 3 are held fixed. 



of a nearby oscillatory state, there was a failure of convergence of the finite-difference code used 
to generate an initial guess for the tangent vector. FigureUf) examines the scaling of the stability 
boundary by plotting pairs of values (A c , 1/Rm). For A c 1 the variation of 1/Rm with A c is 
approximately linear for each curve, with the Ac-dependence weakening for larger A c . Figure \4$ 
also demonstrates that the critical Rayleigh number for convection Rm c (A c ) decreases slightly 
as the concentration ratio increases, so that chimney convection is stable over a wider range of 
(A, Rm) for larger C. The leading order linear variation for A c <C 1 is consistent with a scaling 

A c (Rm) ~ R c /Rm, (3.1) 

for R c constant, which we can reinterpret in terms of an approximate revised stability condition 



Rm.E = RmAc ~ R c , 



(3.2) 
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in terms of an effective Rayleigh number 

Rm_B . gW-ttOlIo (I) , (3 . 3) 
UK \2 J 

formed using the permeability IIo and the horizontal width of a convection cell 1/2. This modified 
Rayleigh number can be rationalized by noting that the mush-liquid interface is a free bound- 
ary and so the dimensional thickness h is free to evolve in response to the flow. Hence, when 
determining a critical porous-medium Rayleigh number gP(Co — C^ITo/t/VK for convection, h 
can scale either with the vertical thermal diffusion length scale k/V, with the convection-cell 
width 1/2, or else some combination of the two. For large Rayleigh number the mushy-layer 
structure is dominated by the convective flow and thus we expect h to be independent of the 
thermal diffusion scale. This leads to h oc 1/2 for Rm S> 1, yielding the effective Rayleigh num- 
ber (|3.3[) . Hence, for strong convection in cells of confined horizontal width, it is the horizontal 
length scale that provides the dominant restriction when determining allowed modes of con- 
vection. This effect is associated with the competition between convection cells for material to 
transport out of the mushy layer. We expect and observe the scaling (|3.2[) to break down when 
A = IV /2k, = O(l), when the vertical thermal length scale becomes comparable to the chimney 
spacing and plays a more significant role in controlling the stability of convecting states. 



4. Scaling of the optimal state for large Rayleigh number 

We now investigate how the states with optimal potential energy fluxes evolve as the vigour of 
convection increases for large Rayleigh number. We begin by identifying relevant scaling regimes 
in the governing equations to provide context for later discussion of the numerical results. 

4.1. Scaling analysis 

We investigate the scaling of the optimal state for strong convection in the asymptotic limit 
Rm 3> 1, with C > 1 and <S > 1 relevant to our numerical simulations. In particular, consistent 
with our hypothesis that the chimney spacing adjusts to provide an optimal solute flux, we derive 
scalings where A varies with Rm (other scalings are possible if the chimney spacing is specified 
a priori and introduces an additional length scale into the problem). Noting that 6 = O(l) as 
a result of the boundary conditions and < <j> < 1 by definition, we determine the possible 
self-consistent asymptotic scalings of (|2.3[) - (|2.16[) that take the form 

i> = Rm'*, x = Rm~ c X, z = RmT d Z, A = Rm _c A, h = Rm~ d #, (4.1) 

capturing convection that penetrates the full depth of the mushy layer. We shall demonstrate 
below that such scalings will only be valid with the further condition Rm/C = O(l), in order to 
maintain a solid fraction with < <f> < 1. With the expectation that X, Z, A, H and * will be 
of order one, the exponents b, c and al can be determined as follows. The balance of conduction 
of heat into the chimney with heat advected along the chimney in (|2.13[) requires c = d + b, 
yielding 

<- Y = »>- (4 - 2 > 

At the mush-liquid interface, the balance of heat advected into the boundary layer with heat 
conducted into the mush requires d — c + b, so that we must have c = d, b — 0, and hence (|2.11[1 
reduces to 

n.V9 = 9«,(v.n-U.n)+of^J, (Z = H), (4.3) 

where V = (d/dX,d/dZ) and U = {-d^/dZ, dV/dX). Next, noting that n = O(l) near the 
mush-liquid interface, a balance of baroclinic torque with viscous dissipation of vorticity in (|2.9|l 
requires c = 1, yielding 



v * - ~ n 8x -( — n— I ■ ^ 
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These conditions result in the asymptotic scalings 

A = °(lii)' h = °{lL)' ^ = °^' for Rm » 1 ' ( 4 - 5 ) 

from which we expect that u = (—dtf)/dz, dip/dx) — O(Rm) and from 1)2.181) the solute flux will 
have the scaling 

F s = O(Rm), (4.6) 

consistent with the observation of I Wells et al\ ()2010l ). 

The rescaled versions of the heat equation (|2.3p and solute equation (|2.4|) become 



^U-W. (4.8) 



JL 
dz 

If U ■ X78 somewhere in the domain, then for Rm ^> 1 (|4.8[) can only be satisfied in a 
self-consistent fashion if Rm/C = O(l), so that C > 1. Hence (|4.8[) predicts a solid fraction scale 

/Rm 



^ = 0^ — )■ (4-9) 

This balance breaks down for Rm 3> C. In particular, we expect that as <j> approaches 1, the per- 
meability hi — > and reduces the strength of flow given by the rescaled vorticity equation (|4.4[) • 
This form of scaling analysis cannot lead to scaling laws suitable for Rm S> C, and it is possible 
that the variation in permeability may lead to convective flow that does not penetrate the full 
depth of the mushy layer. The scaling (|4.9[) can be used to show that (|4.7[) yields a self-consistent 
balance provided S/C = O(l). 

The boundary conditions (I2.15l) - (|2.16|) at z = and x — X retain the same form after rescaling, 
and at the chimney wall (|2.14|l yields 

v '^ e=o {iL)' ( x =°)- ( 4 - iq ) 



and the mass flux condition yields 



Rm 3 
Da 



3 (a i\ 1 9 * 

20 {e + 1) + mdx 



+ aRm|J, (X = 0). (4.11) 



Note that there is upflow with 9* / dX > in the neighbourhood of the chimney, and hence (|4.11j) 
provides two possible consistent scalings for the chimney width valid for two different ranges of 
Rm 2 Da. If Rm 2 Da < 1, then 



.Rm 1 ' 3 , 

and the final term in (|4.11[) is negligible. Poiseuille flow in the chimney is driven by buoyancy 
in the chimney and the pressure gradient exerted by the flow in the mush, both forces being 
of the same asymptotic magnitude. In this limit, the dimensional chimney width 2a oc If = 

[vK,/g/3(Co — Ce)] 1//3 has a length scale If relevant for convection in a purely liquid region. The 
second possible scaling for a occurs when Rm 2 Da S> 1, and we find 

a = o(^), (4.13) 

so that the final term dominates the right-hand side of (|4.11[) so that there is plug flow in 
the chimney. Here the chimney forms a passive conduit of enhanced permeability, adjusting to 
accommodate the flux supplied by convection in the mushy region. The dimensional chimney 
width satisfies 2a oc l p = UK,/g/3(Co — Ce)^Iq consistent with a lengthscale l p controlled by 
porous medium convection in the interior of the mushy layer. The group Rm 2 Da = (lf/l P ) 3 
can therefore be interpreted as a measure of the relative strengths of porous medium convection 
compared to pure fluid convection. We note here that whilst solutions in the limit (|4.13[) remain 
mathematically consistent solutions of mushy-layer theory, the manifestation of the chimney in 




Figure 5. (a) Variation of the optimal solute flux F with mush Rayleigh number Rm, for a 
variety of concentration ratios C. Line styles and symbols are indicated in the legend. All other 
parameters S = 5, (9oo = 0.4 and Da = 5 x 10 -3 are held fixed. For moderate Rm these curves 
can be approximated by F ~ a(Rm — Rmo), where a and Rmo are independent of Rm at 
leading order. (6) Variation of slope a(C) = dF /dKm for three different values of Rm indicated 
in the legend, corresponding to the near- linear sections of the curves in panel (a). Data points 
indicated by symbols are joined by linear interpolation to clarify the underlying trend. 



particular physical contexts may depend on the details of the microstructure of the material in 
question, as discussed in more detail at the end of i |4.2l 

It is interesting to note that if a is determined consistent with the condition (|4.11|) . the 
asymptotic scalings of the remainder of the flow remain independent of the particular scal- 
ing of a. This emphasizes the role of chimneys as passive conduits that assist the drainage of 
fluid driven by convection in the interior of the mushy layer. The scalings (|4.5[l - (|4,6[) . (|4.9p 
and (|4.12|l have also been recovered independently in a recently developed simplified model 
of chimney co nvection that assumes const ant permeability and a linear background tempera- 
ture gradient (jRees Jones fc Worsteril2013D . and we now compare the scalings to full numerical 
simulations where the assumptions of the simplified model are relaxed. 

4.2. Numerical results in the optimal state 

States with optimal fluxes were identified in the numerical simulations for a range of values 
of Rm and C. Guided by the scaling analysis in t j4.ll we now investigate the behaviour of the 
optimal state as Rm varies with all other parameters held fixed. Such variation in the vigour 
of convection could be achieved in directional solidification experiments by varying the growth 
rate V. In this section we fix the par ameter values S = 5, flop — 0.4 and Da = 5 x 10~ 3 for 
consistency with the previous study of IChung fc Wor stcr (2002). 

Figure [5]a shows the variation of the optimal solute flux F with mush- Rayleigh number Rm 
for a wide range of concentration ratios 2 < C < 75. For the largest concentration ratios, the flux 
varies approximately linearly with Rayleigh number across the full plotted range. For smaller 
values of C, initially the flux varies approximately linearly with Rayleigh number before the 
growth begins to saturate with sub-linear growth for larger values of Rm. This is consistent 
with the scaling behaviour identified in i)4.1l with F ~ oRm for Rm 3> 1 and Rm/C = O(l), 
where the slope a is constant, but the scaling breaks down when Rm > C as the solid fraction 
becomes large and reduces the permeability. The magnitude of the solute flux increases with C, 
and in figure 0> we illustrate the C-dependence of the slope a — dF /dKm for three moderate 
values of Rm in the region of near- linear scaling where F ~ a(Rm — Rmo). The increase in 
flux with C is also consistent with variation of the permeability, with the scaling (|4.9[) implying 
smaller solid fraction and hence higher permeability for larger values of C. Figure [6]a provides 
further support for these mechanisms, where we notice that the maximum solid fraction <f> ma x is 
larger for smaller values of C, increases with Rm, and begins to saturate for large values of Rm 
and small C. To further investigate the scaling (|4.9[) . figure [6t> illustrates the variation of the 
maximum solid fraction with Rm/C Data for 4 < C < 75 collapse onto a single curve, with <f> ma x 
varying approximately linearly with Rm/C for small Rm/C, consistent with the scaling (|4.9[) for 
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Figure 6. (a) Variation of the maximum solid fraction (j>max in the domain with mush Rayleigh 
number Rm, for a variety of concentration ratios C. (b) Collapse of the same data plotted as a 
function of Rm/C. Curves for C = 3 and C — 2 deviate slightly from the collapsed trend, which 
is consistent with the hypothesized breakdown of the scalings when C = 0(1). Line styles and 
symbols are indicated in the legend, and are identical to figure [5] All other parameters <S = 5, 
Goo = 0.4 and Da = 5 x 10~ 3 are held fixed. Kinks arising in the curves for small Rm are 
associated with a discontinuous change in spatial position of the maximum value of <f>. 



Rm/C = O(l) and Rm, C 3> 1. Data for C = 2 and C = 3 follow the same qualitative trend 
whilst showing a slight deviation from the collapse of the remaining data. This is consistent with 
the hypothesized breakdown of the scalings when C = O(l). 

Figure [7] demonstrates the scaling behaviour of the optimal chimney half-spacing A and 
corresponding mean mush thickness h, maximum vertical velocity w max and maximum stream- 
function ipmax- Motivated by the scalings (|4.5[) we plot 1/A and \/h versus Rm in panels a 
and c respectively. Both plots show an initial period of linear variation, consistent with the 
scalings A = 0(Rm~ ) and h = 0(Rm _1 ), before the previously discussed flow saturation ef- 
fect leads to deviation from the linear trend at larger Rm. For a given value of Rm both the 
optimal chimney spacing A D and mean mush depth h decrease as C increases. Note that the 
sc alings A oc Rm -1 an d h oc Rm -1 give rise t o the 0(l)-aspect-ratio convective cells identified 
by IWells et al\ (|2010T l and IWells et al\ (1201 if) . and also observed in enthalpy method simula- 
tions in wide cells bv lKatz fc Worsted (|2008T ). The absolute maximum vertical velocity |to m aa;| 
initially increases approximately linearly with Rm (panel b) with stronger flow for larger C. The 
maximum streamfunction ip ma x remains O(l) throughout, with ipmax increasing slightly with 
C. Note that for 25 < C < 75, ipmax increases slightly with Rm after an initial transient, whilst 
for C < 15 we observe that ipmax begins to decrease for large Rm consistent with saturation of 
the mass flux for Rm S> C. 

The behaviour of F 3 , cpmax, A , h, w ma x, and ipmax are consistent with the hypothesized 
scaling regime, which suggests the following picture for the development of the flow in directional 
solidification. As the strength of convection increases, enhanced advection of heat from the 
overlying liquid gives rise to a thinner mushy layer with h oc Rm -1 . Optimal drainage of potential 
energy is achieved with 0(l)-aspect-ratio convection cells which minimize the resistance to flow, 
and hence the optimal chimney spacing also reduces with A oc Rm -1 . The mass flux through 
each convective cell remains of similar magnitude, with ipmax = 0(1) as Rm increases, but the 
smaller chimney spacing allows a higher density of convective cells so that the fluid velocity 
u = O(Rm) and solute flux F 3 = O(Rm) both increase with Rm. The velocities and fluxes 
are smaller for smaller C because the solid fraction increases as C decreases, and the decreased 
permeability provides more resistance to flow. For sufficiently strong convection with Rm S> C 
there is a saturation effect as the solid fraction becomes large, the permeability is reduced, 
and the flow and solute fluxes grow more slowly. It is important to note that this picture for 
directional solidification does not translate directly to the different setting of transient growth, 
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Figure 7. Variation with Rayleigh number Rm of flow properties in the optimal-flux state, 
including (a) inverse-chimney spacing 1/A , (b) absolute maximum flow velocity |io m(M |, (c) 
inverse-thickness of the mushy-layer 1/h, (d) maximum streamfunction value tp m ax as a measure 
of fluid flux. Each curve corresponds to a different value of the concentration ratio C, and the line 
styles and symbols are identical to figure [5] with C = 2 (blue dashed curve with o) , C = 3 (green 
dashed curve with v) > C = 4 (red dashed curve with *), C — 5 (cyan solid curve with >), C = 10 
(magenta solid curve with A), C = 15 (black solid curve with *), C = 25 (blue solid curve with 
o), C — 50 (green solid curve with □), and C = 75 (red solid curve with x). All other parameters 
S = 5, (9oo = 0.4 and Da = 5 x 10~ 3 are held fixed. Kinks arising in the curves for small Rm 
are associated with a discontinuous change in spatial position of the relevant maximum value. 



where h(t) grows over time and the observed coarsening of chi mney spacing as h(t) increases is 
consistent with the aspect ratio X/h = O(l) ({Wells et aUl2010T ). 

In i|4.1l we identified two different possible scalings for the chimney width given by 1)4.121) 
and (|4.13l) . To investigate which of these scalings is relevant we plot a -3 versus Rm in figure [HJa 
and a" 1 versus Rm in figure [Sf>. For large values of Rm, figure [S] shows that a" 3 increases 
nonlinearly with Rm, whilst a -1 shows approximately linear variation with Rm. This suggests 
that for large Rm our numerical simulations with Da = 5 x 10 -3 correspond to the scaling 
a — <9(1/Rm), consistent with the limit Rm 2 Da 3> 1. However, one would expect the alternative 
scaling (|4.12|) to be achieved for smaller values of Da. For the simulations with 2 < C < 5, both 
1/a 3 and 1/a diverge for small Rm, with the chimney closing and a — > as the Rayleigh number 
approaches the critical threshold for existence of chimney convection. 

We conclude this section by returning to the criterion for physical applicability of ideal mushy 
layer theory. The scaling behaviours and numerical solutions discussed above are mathemati- 
cally consistent solutions of the relevant ideal mushy layer equations, which contain no direct 
information about the material microstructure. However, the details of this microstructure may 
affect the applicability of certain aspects of the theory to particular physical settings. We focus 
our discussion on materials where the permeability n o ~ S 2 scales as the square of the dimen- 
sional pore scale 5 (neglecting any pre-factor in the scaling), so that the dimensionless pore 
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Figure 8. Comparison of different scaling regimes for the chimney width a in the optimal state 
as a function of Rayleigh number Rm, with each curve representing a different concentration ratio 
C indicated by different line and symbol styles. For large Rm, panel (a) shows that 1/a 3 increases 
nonlinearly with Rm, whilst panel (6) indicates that 1/a scales approximately linearly with Rm, 
consistent with the scaling (|4.13[) . Curves for small Rm and 2 < C < 5 show a divergence of 
1/a 3 and 1/a, because a — > and the chimneys close as the Rayleigh number approaches the 
critical threshold for existence of chimney convection. The line styles and symbols are identical 
to figure [5] with C — 2 (blue dashed curve with o) , C = 3 (green dashed curve with v) > C — 4 
(red dashed curve with *), C — 5 (cyan solid curve with >), C = 10 (magenta solid curve with 
A), C = 15 (black solid curve with *), C — 25 (blue solid curve with o), C = 50 (green solid 
curve with □), and C = 75 (red solid curve with x). All other parameters 5 = 5, (9oo = 0.4 and 
Da = 5 x 10~ :! are held fixed. 



scale satisfies S ~ Da 1 / 2 . For this class of materials, we can estimate the pore-space Reynolds 
number as Re p ~ US/v ~ \w\ Da 1//2 /cr where a = u/k is the Prandtl number. The asymp- 
totic scalings (|4.5p yield Re p = O ^RmDa 1//2 /<jj , so that the Darcy flow approximation breaks 

down when RmDa 1 ' 2 ^> a. Note that the parameter range of our numerical results leads to 
\w\ Da 1//2 < 0.5, so that the Darcy flow approximation is reasonable for any fluids of moderate 
to large Prandtl number. The physical details of the microstructure may also be significant for 
determining whether macroscopic flow features are larger than the pore scale. The chimneys 
provide the smallest macroscopic lengthscale in the system, with two different scalings. The 

scaling (|4.12p results in a ratio S/a — O LRm^Da 1 / 8 ^ so that the chimney width is wider 

than the pore scale in the asymptotic limit RmDa 1 ' 2 — > 0. By contrast, the scaling 1)4.131) 

gives S/a = O ^RmDa 1 '^ so that materials with n o ~ S 2 have a predicted chimney width 

smaller than the physical pore scale for sufficiently large values of RmDa 1 ' 2 . Hence, whilst the 
mathematical description remains self consistent for RmDa 1//2 2> 1, the physical manifestation 
of the chimney shape in the limit (|4. 13|) may depend on the permeability and microstructure 
of the material in question. However, we emphasize that the asymptotic flow structure in the 
interior of the mushy layer remains relatively insensitive to the particular shape of the chimney 
with the scalings (|4.5|) . (|4.6[) . and (|4.9[) for the interior of the mushy layer remaining valid in 
both regimes (|4.12p and (|4.13|) for the chimney. In particular, the value of Da only enters the 
problem through ()2.12[) , and hence if chimney scalings are determined consistent with (|2.12l) we 
expect the dominant structure of the flow in the interior of the mushy layer to be robust, and 
independent of both the particular chimney shape and value of Da . This would be the case if 
one were to solve the heat equation (|2.3p and the vorticity equation 1)2.91) subject to boundary 
conditions (I2.13|) and (|2.14[) . with a adjusted to satisfy (|2.12|) . The physics of the initial problem 
and our numerical methodology correspond to a solution of the heat equation (|2.3p and vor- 
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Figure 9. Comparison of stable steady states of solute flux F s as a function of chimney half-s- 
pacing A using (a) the cubic permeability function (|2,5|l . and (b) the Carman-Kozeny (C-K) 
permeability function (|2,6|l using e = 0.1. Both permeability functions result in hysteretic be- 
haviour, with an upper stable branch of chimney convection for A > A c (black curve) and a lower 
branch with a state of no flow being stable for A < \ u (blue line with triangles) . The parameters 
Rm = 35, C = 15, S — 3, 8^, = 0.5 and Da = 10 -4 are held fixed in both plots. 

ticity equation (|2.9[) using boundary conditions (|2.13|l and (|2.12[) at x = 0, with a adjusted to 
satisfy the marginal equilibrium condition (|2.14|l . Whilst mathematically rational, we have not 
proved that these two different approaches are guaranteed to give the same solution. However, 
it is interesting to note that previous simulations that varied Da with all other par ameters held 
fixed resulted in a mass flux through the chimney wall that was independent of Da ()Zhong et al\ 
I2013) . 



5. Dependence of flow on the permeability function 

We conclude the presentation of results by discussing a selection of simulations using the 
Carman-Kozeny permeability function (|2.6[) to examine the influence of the permeability varia- 
tion on the flow dynamics. The stable solution branches of F S (X) using the cubic permeability 
function (|2.5[) are compared in figure [9]a. with those using the Carman-Kozeny permeability 
function (|2.6|l in figure \§}>, with all other common parameters the same in both plots. The two 
different permeability functions produce the same qualitative structure of flow dynamics, with 
the usual hysteresis loop involving a state of no flow for A < \ u (blue line with triangles) and 
an upper branch of chimneys convection for A > A c (black solid curve). This indicates that 
the overall qualitative structure of the flow dynamics is robust for both of the different perme- 
ability functions. The main quantitative differences are that the Carman-Kozeny permeability 
permits a stronger flow with larger solute fluxes, and the bifurcation points A c and X u are at 
smaller A which indicates that chimney convection is more favourable for the Carman-Kozeny 
permeability than for the cubic permeability. 

The optimal solute flux also occurs at a smaller chimney half-spacing A for the Carman- 
Kozeny permeability and in figure \W\ we compare the corresponding flow profiles for cubic 
permeability (panels a-c) and Carman-Kozeny permeability (panels d-f). The solid fraction <j> 
shows greater variation for the Carman-Kozeny permeability than the cubic permeability (see 
colour scales in figures ITul s and b respectively), which will result in a greater contrast in perme- 
ability Tl(4>) across the mushy region. This suggests flow focussing as a possible mechanism to 
explain the difference in optimal chimney spacing, with the larger permeability contrast for the 
Carman-Kozeny flow encouraging the flow to concentrate in narrower convection cells in order 
to reduce the resistance to flow. The narrower chimney spacing generates a larger baroclinic 
torque RmII<9# / dx = 0(Rm/A) in the vorticity equation (|2.9[) which drives stronger flow. This 
is consistent with flow streamlines for the Darcy velocity u in figures ITOla.d and net fluid flux 
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Figure 10. Comparison of steady-state mushy layer profiles with optimal chimney spacing 
A = A for cubic permeability function (|2.5[) (panels a-c) and Carman-Kozeny permeability 
function (|2.6p (panels d-f). The colour scales show dimensionless temperature 6(x,z) in panels 
(a, d), and solid fraction (f>(x,z) in panels (6, e) with black solid contours of constant 9 and 
<j> respectively. Streamlines of Darcy velocity u are shown by magenta curves with increment 
A 1 !/ = 0.25 in panels (a,d), and streamlines of the net fluid flux q = u k are illustrated in 
panels (6,e). In addition, the dashed magenta curves in panels (a,b) show intermediate contours 
with increment A\l/ = 0.125. Panels (c, /) illustrate the variation of the chimney width a(z) 
with height. All other parameters are identical to figure [2] 

q = u — k in figures 110b . e, both of which indicate stronger flow for the Carman-Kozeny per- 
meability. The stronger flow enhances vertical heat transport through the chimney and results 
in the isotherms having steeper slope at the chimney wall consistent with the boundary condi- 
tion (|2.13|) . This results in the larger curvature of isotherms observed for the Carman-Kozeny 
permeability ffigurefTObO than for the cubic permeability (figure [TObV The larger isotherm slope 
for the Carman-Kozeny permeability also contributes to increasing the baroclinic torque and 
driving stronger flow. Note that the mush-liquid interface is an isotherm, and also shows sig- 
nificant curvature for the Carman-Kozeny permeabili ty in a manner remin iscent of previous 
enthalpy method simulations (compare to figure 10 of iKatz fc W orstcr 2008). The chimney is 
wider for the Carman-Kozeny permeability in order to accommodate the larger flux of fluid 
escaping the mushy layer. 



6. Conclusions 

We have considered finite-amplitude convection in a mushy layer with a periodic array of 
chimneys formed during directional solidification, and studied the parametric dependence of the 
flow over a wide range of chimney spacings, Rayleigh numbers and concentration ratios. A selec- 
tion of simulations with a Carman-Kozeny permeability ()2.6|) illustrate that the flow dynamics is 
qualitatively similar to flow with a cubic permeability function (|2.5p . with the Carman-Kozeny 
permeability resulting in stronger flow and narrower chimney spacings, potentially consistent 
with the process of flow focussing. For the cubic permeability function and a fixed concentra- 
tion ratio, stable states of chimney convection exist for sufficiently large chimney spacings and 
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sufficiently large Rayleigh number, consistent with a stability boundary (|3 . 2p characterised by 
the effective Rayleigh number 1)3. 3[) based on the convection cell width 1/2 when IV/2k <C 1. 
This recognizes the role of the horizontal wavelength in controlling mode restriction for cells 
that are narrower than the vertical thermal length scale. As the concentration ratio increases, 
the observed states of chimney convection are s table over a wider range of A and Rm. 

The optimal potential energy flux criterion of lWells et al\ (|2010T ) was hypothesized as a phys- 
ical mechanism to select preferred chimney spacings observed during directional solidification, 
and the scaling behaviour of these optimal flux states has been investigated here. For moder- 
ately large Rayleigh numbers Rm ^> 1 with Rm/C = 0(1), scaling laws (|4.5p . (|4.6p and (|4.9jl 
have been developed consistent with numerical simulations. These scalings illustrate that as the 
strength of convection increases, the flow velocity increases by generating a larger number of 
narrower convection cells of O(l) aspect ratio in a thinner mushy region. The scalings derived 
here are also observed in a recent simpli fied model that assumes con stant permeability and a 
linear background temperature gradient ()Rees Jones fc Wor stcr 2013|), which provides support 
that their approximations do not have a major influence on the flow scaling. The flow dynamics 
in the chimney exhibits two asymptotic regimes: buoyancy and pressure gradients imposed from 
the neighbouring mush are of similar magnitude for Rm 2 Da <C 1 producing the scaling (|4.12p . 
whilst for Rm 2 Da 2> 1 the chimney forms an entirely passive conduit that adjusts in width to 
accommodate the fluid flux from the mushy region, giving the scaling (|4.13[1 , For Rm/C = O(l) 
the scaling laws rec over a solute flux t h at increases l inearl y with Rayleigh number, consistent 
with the results of IWells et al\ (|2010f ). IWells et al\ (|201ir i. and a recent approximate model 
l)Rees Jones fc W orstcr 20131). However, the linear scaling does not persist indefinitely and for 
Rm > C we observe sub-linear growth resulting from saturation of the flow due to porosity- 
dependent permeability as the solid fraction increases. This flow saturation for Rm 2> C raises 
interesting questions regarding the applications in sea ice modelling, where typically C C 1. 
Intriguingly, IWells et al\ (|20 1 ll ) found that the salinity evolution in experimental sea ice growth 
was consistent with a scaling F s ~ a(Rm — Rmo), despite being in the regime C <C 1. It remains 
an open question to determine whether directional solidification will again recover F s cx Rm in 
a new dynamical regime for C C 1, or whether this is associated with some specific property 
of transient solidification where the ice thickness h increases over time. A further question of 
interest relates to how convection in the neighbouring liquid region, which was neglected in this 
study, interacts with the mushy layer flow dynamics described here. If there is turbulent convec- 
tion in the liquid, one might intuitively expect a thermal and compositional boundary layer to 
develop ahead of the mush-liquid interface, so that the interfacial temperature and concentra- 
tion differ from those of the far-field liquid. One would also need to account for modified heat 
and solute fluxes at the mush-liquid interface, and the precise details of this interaction remain 
to be explored. Nevertheless, it is hoped that the scaling behaviour and system properties illus- 
trated here may provide insight for simplified modelling applicable to both sea ice growth and 
other problems involving mushy-layer convection. 
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Appendix A. Heat flux at the chimney wall 

In this appendix we show that the boundary condition (|2.13|l for the heat flux at the chimney 
wall, previously derived by ISchulze fc Worsted (|l998f ) in the asymptotic limit Rm 4 ^ 3 Da <C 1, 
can be justified as a perturbation approximation under the ansatz a/h <C 1 and ip = O(l) for a 
wider range of Rm and Da . Motivated by this ansatz, we use rescaled coordinates £ = x/eh(0, t), 
£ = z/h(0,t) to resolve flow within the chimney, and a rescaled chimney width A = a/eh(0,t). 
It is assumed that £, £ and A are all of order one and e < 1, so that the chimney occupies 
< £ < A and a/h <C 1. If tp — 0(1), then in rescaled co-ordinates the heat equation (|2,3[) 
within the chimney reduces to the leading-order balance 

= O(e). (Al) 

Integrating twice and applying the symmetry condition 80/d£, = at £ = yields — f{C)+0(e) 
for some function /(C), so that 6 is independent of £ at leading order. 
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To determine the heat flux at the chimney wall, we consider the first order correction to (|A 
with (|2.3[) yielding 

d 2 9 dtp 86 2 , 

W = e idt + oie) ' (A2) 

where we have taken advantage of the fact that d8/d£ = O(e). Noting that 89 /dC, is independent 
of £ at leading order, integrating (IA 2 [I once and using the boundary conditions d6/d^ = -0 = 
at £ = yields the condition 

|=^| + 0(e 2 ), (A3) 



which is valid for < £ < A- Evaluating (|A 3|l &t ^ = A, and then returning to (a;, z) co-ordinates 
recovers the heat flux condition (|2.13[) at the chimney wall. 



Appendix B. Notes on the numerical procedure 

Two complementary numerical methods were applied, as described in the main text. We 
elaborate on the time-dependent code in tffi.ll and the arc-length continuation scheme in . 2 1 

B.l. Time- dependent solution 

The time-dependent problem was solved using finite differencing. We use the co-ordinate trans- 
form £ = z/h to allow the use of a time-independent spatial grid with < £ < 1. Semi-implicit 
Crank-Nicholson timestepping is used to update the enthalpy 9 — <f>S /C and the total concentra- 
tion C(f> + 8(1 — <j>), using streamfunction data from the previous timeste p only. The tridiagona l 
system resulting from (|2.4p was solved using standard LAPACK routines (| Anderson et aZ.lll999f ). 
After time discretization, (|2.3|) yields an elliptic system in space which, along with the vorticity 
equation (|2.9[), was solved u sing multigrid iteration via an adaptation of the MUDPACK multi- 
grid routines IjAd ams 198^). The boundary condition (|2.12[) for the vorticity equation for the 
updated streamfunction ip" +1 is split in the form 

V>" +1 = (1 - ei At)ip n + eiAt 

where £i is constant and At the numerical time step. This relaxation approach satisfies (|2.12|l 
with a discretization error of O(At). Using streamfunction data ip n from the previous timestep 
allows a straightforward implementation of the MUDPACK multigrid routines when a = over 
a subsection of the boundary. The modified free-boundary condition (|2.11[) and (|2.17[) is solved 
by minimizing the least-squares residual from the relation 

[ {n-Vfl + fooKu-kj-n-V-n]} 2 dx = 0, (B2) 
oho Jo 

which gives an algebraic equation for a new estimate h = h* at each timestep. The updated 
values of h a and a(z,t) are then calculated using relaxation with 

The choice of the relaxation rates e\ , €2, and £3 used here influence the solution behaviour over 
very short timescales, but do not influence the final steady state except for close approach to the 
bifurcation points, where arc-length continuation was used to confirm the solution trajectory. 
Values ei < 30, £2 < 20, and £3 < 5 were found to give good numerical convergence, with the 
maximal values typically used to obtain rapid relaxation compared to the development of the 
flow. Each subsection of the computer code was tested independently for accurate convergence 
against analytic solutions of simple advection and diffusion problems, and the full code was 
tested by adding false forcing to generate known analytic solutions of the modified system. The 
final solutions were tested by reducing step size Aa;, and checking for 0(Ax 2 ) convergence. 

B.2. Arc-length continuation 

A pseudo-arcle ngth continu ation method was used to trace the branches of stable and unstable 
steady states ijKellerl 119771 ). Using N x x spatial grid points, the finite-difference solution 



3Da(l - 
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has on the order of N = 3iVj - + 2 de grees of freedom accounting for the variables 0(x,z), 
<f>(x, z), ijj(x, z) and a(z), along with ho and hi which characterize h(x, t), and subtracting known 
9, <f> and tjj values specified at boundaries. Rather than evaluating and performing inversions for 
the resulting large dense N x N Jacobian matrix, in an analogue with pseudo-spectral methods 
we seek a lower order representation of the solution expanded in Chebyshev polynomials 

M,M M 

(<9,0,VO= (0™«>^Vw)T m (2|-l)T„(2C-l), a{z) = ^a n T n (2C- 1), (B4) 

n,m=0 n — 

where the Chebyshev polynomials are conveniently defined by 

T„(cos 6») = cos n9 for < 9 < tt. (B 5) 



As a result of the identity (|B 5[) . the finite-difference solution is straightforwardly projected onto 
the Chebyshev polynomials using Fast Fourier Transforms. 

The pseudo-arclength continuation scheme, which is an application of the method of I Keller! 
(| 19771 ). proceeds as follows. We symbolically represent the finite difference version of (I2.3I) - (|2.16[) 

as 

§T=/(y,7), (B6) 

where y = {9, cj>, ip, a, ho, hi} represents the finite-difference solution and 7 is the parameter to 
be varied. We let V denote projections into Chebyshev space according to (|B 4[) with V^ 1 the 
corresponding inversion. Two initial steady states yo and y_i, with /(yo, 70) = /(y-i, 7-i) = 0, 
form an initial approximation to the tangent vector in (Vy, 7) space viz., 

t (Pyo-Py-1,70-7-1) mr , 
° ||(Pyo-Py-i,7o-7-i)ir 1 ' 

We seek a new solution with pseudo-arclength step As by using the predictor step (Vyijji) = 
(Vyo, 70) + Asto, and a corrector based on Newton iteration to solve the coupled system 

?V(yo,7o) = 0, (B8) 

to ■ (Ty - Vyo, 7 - To) = As, (B 9) 

for the new solution (Vy, 7) in Chebyshev space. Note that all the necessary function evaluations 
of Vf(y, 7) are carried out sequentially as follows. First, we transfer the Chebyshev projection 
Vy back to finite difference space using y = V~ Vy. Secondly, we compute the time evolution 
of (|2,3|1 - (|2,16[) with a(z) held fixed to estimate f(y, 7). Finally, we return to the Chebyshev space 
via Vf(y, 7) with the residual in satisfying (|2.14p used in place of estimates of da/dt. After the 
corrector step yields a new solution estimate (Vy,^), a steady version of the finite-difference 
code is used to post-process the solution. Iteratively repeating the above predictor-corrector 
process allows the trajectory of steady-state solution branches to be traced. 
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